Australia LCA

Below shows:

  1. Best model item response probability
  2. 3D plot of class-conditional item response probs

min_bic <- 100000
for(i in 2:7){
  lc <- poLCA(f, AUS, nclass=i, maxiter=3000, 
              tol=1e-5, na.rm=FALSE,  
              nrep=10, verbose=TRUE, calc.se=TRUE)
  if(lc$bic < min_bic){
    min_bic <- lc$bic
    LCA_best_model<-lc
  }
} 
print(LCA_best_model)
Conditional item response (column) probabilities,
 by outcome variable, for each class (row) 
 
$tax
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.1756 0.0208 0.0903 0.0288 0.1086 0.0753 0.0965 0.0741 0.0437 0.2863
class 2:  0.0257 0.0288 0.0595 0.1047 0.3553 0.1350 0.1489 0.0897 0.0346 0.0178
class 3:  0.0028 0.0303 0.0672 0.0739 0.1374 0.1453 0.2182 0.1862 0.1188 0.0199

$religion
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.7301 0.0227 0.0434 0.0198 0.0756 0.0169 0.0080 0.0205 0.0110 0.0521
class 2:  0.1514 0.1858 0.1904 0.0893 0.2203 0.0724 0.0563 0.0152 0.0144 0.0046
class 3:  0.4164 0.1947 0.1566 0.0888 0.0643 0.0328 0.0255 0.0169 0.0042 0.0000

$free_election
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.0496 0.0024 0.0020 0.0000 0.0221 0.0060 0.0000 0.0059 0.0112 0.9007
class 2:  0.0208 0.0360 0.0430 0.0392 0.1787 0.0830 0.1968 0.1702 0.1254 0.1068
class 3:  0.0018 0.0000 0.0025 0.0008 0.0000 0.0063 0.0186 0.0880 0.2217 0.6604

$state_aid
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.1045 0.0245 0.0472 0.0400 0.1195 0.0384 0.0310 0.1055 0.0360 0.4533
class 2:  0.0497 0.0778 0.0712 0.0698 0.3181 0.1734 0.1217 0.1004 0.0168 0.0011
class 3:  0.0223 0.0355 0.0650 0.0691 0.0984 0.1346 0.2168 0.1794 0.1449 0.0340

$civil_rights
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.0921 0.0051 0.0100 0.0174 0.0759 0.0152 0.0260 0.0449 0.0229 0.6906
class 2:  0.0294 0.0765 0.0656 0.0827 0.3384 0.1512 0.1644 0.0917 0.0000 0.0000
class 3:  0.0000 0.0068 0.0168 0.0289 0.0399 0.0527 0.0915 0.2350 0.2678 0.2606

$women
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.0316 0.0009 0.0000 0.0000 0.0140 0.0021 0.0036 0.0025 0.0092 0.9361
class 2:  0.0116 0.0109 0.0111 0.0185 0.1493 0.0478 0.1579 0.2057 0.0992 0.2879
class 3:  0.0000 0.0049 0.0000 0.0000 0.0042 0.0000 0.0009 0.0528 0.1769 0.7603

Estimated class population shares 
 0.3989 0.2019 0.3991 
 
Predicted class memberships (by modal posterior prob.) 
 0.3967 0.1961 0.4072 
 
========================================================= 
Fit for 3 latent classes: 
========================================================= 
number of observations: 1336 
number of estimated parameters: 164 
residual degrees of freedom: 1172 
maximum log-likelihood: -12832.31 
 
AIC(3): 25992.62
BIC(3): 26845
G^2(3): 8045.548 (Likelihood ratio/deviance statistic) 
X^2(3): 4059074 (Chi-square goodness of fit) 
 
plot(LCA_best_model)


For the Sweden data

Below shows:

  1. Best model item response probability
  2. 3D plot of class-conditional item response probs
  3. 2D profile plots of mean response
min_bic <- 100000
for(i in 2:7){
  lc <- poLCA(f, SWEDEN, nclass=i, maxiter=3000, 
              tol=1e-5, na.rm=FALSE,  
              nrep=10, verbose=TRUE, calc.se=TRUE)
  if(lc$bic < min_bic){
    min_bic <- lc$bic
    LCA_best_model<- lc
  }
} 
LCA_best_model
Conditional item response (column) probabilities,
 by outcome variable, for each class (row) 
 
$tax
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:  0.0134 0.0553 0.1039 0.0425 0.1645 0.1149 0.2246 0.2117 0.0599 0.0093
class 2:  0.0674 0.0549 0.0691 0.0660 0.0739 0.0572 0.1369 0.1936 0.1058 0.1751

$religion
           Pr(1)  Pr(2)  Pr(3)  Pr(4)
class 1:  0.3150 0.2601 0.1875 0.2374
class 2:  0.8031 0.0864 0.0456 0.0649

$free_election
          Pr(1) Pr(2) Pr(3) Pr(4) Pr(5) Pr(6) Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:      0     0     0     0     0     0     0 0.3638 0.3458 0.2904
class 2:      0     0     0     0     0     0     0 0.0217 0.0166 0.9617

$state_aid
           Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)  Pr(6)  Pr(7)  Pr(8) Pr(9) Pr(10)
class 1:  0.0125 0.0358 0.0725 0.0491 0.0912 0.1206 0.1678 0.2744 0.120 0.0562
class 2:  0.0704 0.0420 0.0395 0.0362 0.0678 0.0595 0.0866 0.1729 0.126 0.2991

$civil_rights
          Pr(1) Pr(2) Pr(3) Pr(4) Pr(5) Pr(6)  Pr(7)  Pr(8)  Pr(9) Pr(10)
class 1:      0     0     0     0     0     0 0.2511 0.3321 0.3198 0.0970
class 2:      0     0     0     0     0     0 0.0245 0.0460 0.0722 0.8573

$women
          Pr(1) Pr(2) Pr(3) Pr(4) Pr(5) Pr(6) Pr(7) Pr(8) Pr(9) Pr(10)
class 1:      0     0     0     0     0     0     0     0 0.441  0.559
class 2:      0     0     0     0     0     0     0     0 0.061  0.939

Estimated class population shares 
 0.249 0.751 
 
Predicted class memberships (by modal posterior prob.) 
 0.2417 0.7583 
 
========================================================= 
Fit for 2 latent classes: 
========================================================= 
number of observations: 964 
number of estimated parameters: 97 
residual degrees of freedom: 867 
maximum log-likelihood: -6750.805 
 
AIC(2): 13695.61
BIC(2): 14168.11
G^2(2): 2635.129 (Likelihood ratio/deviance statistic) 
X^2(2): 14270.99 (Chi-square goodness of fit) 
 
plot(LCA_best_model)

Profile plot for Sweden data (2-classes)

plot = ggplot(profile_long, aes(x = variable, y = value, group = class, color = class)) +
  geom_point(size = 2.25)+
  geom_line(size = 1.25) +
  labs(x = NULL, y = "Mean value of the response", main = "Profile plot") +
  theme_bw(base_size = 14)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "top")
p = ggplotly(plot, tooltip = "all") %>%
  layout(legend = list(orientation = "h", y = 1.2))
print(p)
NULL
LS0tCnRpdGxlOiAiTENBX3ByZWxpbWluYXJ5IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgoKIyNBdXN0cmFsaWEgTENBCkJlbG93IHNob3dzOiAKCjEuIEJlc3QgbW9kZWwgaXRlbSByZXNwb25zZSBwcm9iYWJpbGl0eQoyLiAzRCBwbG90IG9mIGNsYXNzLWNvbmRpdGlvbmFsIGl0ZW0gcmVzcG9uc2UgcHJvYnMKCmBgYHtyfQoKbWluX2JpYyA8LSAxMDAwMDAKZm9yKGkgaW4gMjo3KXsKICBsYyA8LSBwb0xDQShmLCBBVVMsIG5jbGFzcz1pLCBtYXhpdGVyPTMwMDAsIAogICAgICAgICAgICAgIHRvbD0xZS01LCBuYS5ybT1GQUxTRSwgIAogICAgICAgICAgICAgIG5yZXA9MTAsIHZlcmJvc2U9VFJVRSwgY2FsYy5zZT1UUlVFKQogIGlmKGxjJGJpYyA8IG1pbl9iaWMpewogICAgbWluX2JpYyA8LSBsYyRiaWMKICAgIExDQV9iZXN0X21vZGVsPC1sYwogIH0KfSAKCgpgYGAKCmBgYHtyfQpwcmludChMQ0FfYmVzdF9tb2RlbCkKYGBgCmBgYHtyfQpwbG90KExDQV9iZXN0X21vZGVsKQpgYGAKCi0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KIyNGb3IgdGhlIFN3ZWRlbiBkYXRhCkJlbG93IHNob3dzOiAKCjEuIEJlc3QgbW9kZWwgaXRlbSByZXNwb25zZSBwcm9iYWJpbGl0eQoyLiAzRCBwbG90IG9mIGNsYXNzLWNvbmRpdGlvbmFsIGl0ZW0gcmVzcG9uc2UgcHJvYnMKMy4gMkQgcHJvZmlsZSBwbG90cyBvZiBtZWFuIHJlc3BvbnNlCgoKICAKYGBge3J9Cm1pbl9iaWMgPC0gMTAwMDAwCmZvcihpIGluIDI6Nyl7CiAgbGMgPC0gcG9MQ0EoZiwgU1dFREVOLCBuY2xhc3M9aSwgbWF4aXRlcj0zMDAwLCAKICAgICAgICAgICAgICB0b2w9MWUtNSwgbmEucm09RkFMU0UsICAKICAgICAgICAgICAgICBucmVwPTEwLCB2ZXJib3NlPVRSVUUsIGNhbGMuc2U9VFJVRSkKICBpZihsYyRiaWMgPCBtaW5fYmljKXsKICAgIG1pbl9iaWMgPC0gbGMkYmljCiAgICBMQ0FfYmVzdF9tb2RlbDwtIGxjCiAgfQp9IApgYGAKCmBgYHtyfQpMQ0FfYmVzdF9tb2RlbAoKYGBgCgpgYGB7cn0KcGxvdChMQ0FfYmVzdF9tb2RlbCkKYGBgCgpQcm9maWxlIHBsb3QgZm9yIFN3ZWRlbiBkYXRhICgyLWNsYXNzZXMpCmBgYHtyfQoKcGxvdCA9IGdncGxvdChwcm9maWxlX2xvbmcsIGFlcyh4ID0gdmFyaWFibGUsIHkgPSB2YWx1ZSwgZ3JvdXAgPSBjbGFzcywgY29sb3IgPSBjbGFzcykpICsKICBnZW9tX3BvaW50KHNpemUgPSAyLjI1KSsKICBnZW9tX2xpbmUoc2l6ZSA9IDEuMjUpICsKICBsYWJzKHggPSBOVUxMLCB5ID0gIk1lYW4gdmFsdWUgb2YgdGhlIHJlc3BvbnNlIiwgbWFpbiA9ICJQcm9maWxlIHBsb3QiKSArCiAgdGhlbWVfYncoYmFzZV9zaXplID0gMTQpKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gNDUsIGhqdXN0ID0gMSksIGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiKQpwID0gZ2dwbG90bHkocGxvdCwgdG9vbHRpcCA9ICJhbGwiKSAlPiUKICBsYXlvdXQobGVnZW5kID0gbGlzdChvcmllbnRhdGlvbiA9ICJoIiwgeSA9IDEuMikpCnByaW50KHApCgpgYGAKCg==